{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5175177c",
   "metadata": {},
   "source": [
    "# Chapter 1 - Introduction (Python Code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e0179da",
   "metadata": {
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [],
   "source": [
    "! pip install --upgrade quantecon_book_networks kaleido"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36c17bdf",
   "metadata": {},
   "source": [
    "We begin by importing the `quantecon` package as well as some functions and data that have been packaged for release with this text."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf584c62",
   "metadata": {},
   "outputs": [],
   "source": [
    "import quantecon as qe\n",
    "import quantecon_book_networks\n",
    "import quantecon_book_networks.input_output as qbn_io\n",
    "import quantecon_book_networks.data as qbn_data\n",
    "import quantecon_book_networks.plotting as qbn_plot\n",
    "ch1_data = qbn_data.introduction()\n",
    "default_figsize = (6, 4)\n",
    "export_figures = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "858f6c0d",
   "metadata": {},
   "source": [
    "Next we import some common python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "622fff13",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.ticker as ticker\n",
    "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n",
    "import matplotlib.patches as mpatches\n",
    "import plotly.graph_objects as go\n",
    "quantecon_book_networks.config(\"matplotlib\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ee79c32",
   "metadata": {},
   "source": [
    "## Motivation\n",
    "\n",
    "### International trade in crude oil 2021\n",
    "\n",
    "We begin by loading a `NetworkX` directed graph object that represents international trade in crude oil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e495966",
   "metadata": {},
   "outputs": [],
   "source": [
    "DG = ch1_data[\"crude_oil\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "284d4d87",
   "metadata": {},
   "source": [
    "Next we transform the data to prepare it for display as a Sankey diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c944e03",
   "metadata": {},
   "outputs": [],
   "source": [
    "nodeid = {}\n",
    "for ix,nd in enumerate(DG.nodes()):\n",
    "    nodeid[nd] = ix\n",
    "\n",
    "# Links\n",
    "source = []\n",
    "target = []\n",
    "value = []\n",
    "for src,tgt in DG.edges():\n",
    "    source.append(nodeid[src])\n",
    "    target.append(nodeid[tgt])\n",
    "    value.append(DG[src][tgt]['weight'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cd4537f",
   "metadata": {},
   "source": [
    "Finally we produce our plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6a9a856",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = go.Figure(data=[go.Sankey(\n",
    "    node = dict(\n",
    "      pad = 15,\n",
    "      thickness = 20,\n",
    "      line = dict(color = \"black\", width = 0.5),\n",
    "      label = list(nodeid.keys()),\n",
    "      color = \"blue\"\n",
    "    ),\n",
    "    link = dict(\n",
    "      source = source,\n",
    "      target = target,\n",
    "      value = value\n",
    "  ))])\n",
    "\n",
    "\n",
    "fig.update_layout(title_text=\"Crude Oil\", font_size=10, width=600, height=800)\n",
    "if export_figures:\n",
    "    fig.write_image(\"figures/crude_oil_2021.pdf\")\n",
    "fig.show(renderer='svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "663c1374",
   "metadata": {},
   "source": [
    "### International trade in commercial aircraft during 2019\n",
    "\n",
    "For this plot we will use a cleaned dataset from \n",
    "[Harvard, CID Dataverse](https://dataverse.harvard.edu/dataverse/atlas)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31cf9ca7",
   "metadata": {},
   "outputs": [],
   "source": [
    "DG = ch1_data['aircraft_network']\n",
    "pos = ch1_data['aircraft_network_pos']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "990418cb",
   "metadata": {},
   "source": [
    "We begin by calculating some features of our graph using the `NetworkX` and\n",
    "the `quantecon_book_networks` packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4da4ce46",
   "metadata": {},
   "outputs": [],
   "source": [
    "centrality = nx.eigenvector_centrality(DG)\n",
    "node_total_exports = qbn_io.node_total_exports(DG)\n",
    "edge_weights = qbn_io.edge_weights(DG)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ead51019",
   "metadata": {},
   "source": [
    "Now we convert our graph features to plot features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9be3605",
   "metadata": {},
   "outputs": [],
   "source": [
    "node_pos_dict = pos\n",
    "\n",
    "node_sizes = qbn_io.normalise_weights(node_total_exports,10000)\n",
    "edge_widths = qbn_io.normalise_weights(edge_weights,10)\n",
    "\n",
    "node_colors = qbn_io.colorise_weights(list(centrality.values()),color_palette=cm.viridis)\n",
    "node_to_color = dict(zip(DG.nodes,node_colors))\n",
    "edge_colors = []\n",
    "for src,_ in DG.edges:\n",
    "    edge_colors.append(node_to_color[src])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b44e9249",
   "metadata": {},
   "source": [
    "Finally we produce the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dfccd9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "ax.axis('off')\n",
    "\n",
    "nx.draw_networkx_nodes(DG, \n",
    "                        node_pos_dict, \n",
    "                        node_color=node_colors, \n",
    "                        node_size=node_sizes, \n",
    "                        linewidths=2, \n",
    "                        alpha=0.6, \n",
    "                        ax=ax)\n",
    "\n",
    "nx.draw_networkx_labels(DG, \n",
    "                        node_pos_dict,  \n",
    "                        ax=ax)\n",
    "\n",
    "nx.draw_networkx_edges(DG, \n",
    "                        node_pos_dict, \n",
    "                        edge_color=edge_colors, \n",
    "                        width=edge_widths, \n",
    "                        arrows=True, \n",
    "                        arrowsize=20,  \n",
    "                        ax=ax,\n",
    "                        node_size=node_sizes, \n",
    "                        connectionstyle='arc3,rad=0.15')\n",
    "\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/commercial_aircraft_2019_1.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70e419f0",
   "metadata": {},
   "source": [
    "## Spectral Theory\n",
    "\n",
    "### Spectral Radii\n",
    "\n",
    "Here we provide code for computing the spectral radius of a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0ac0f17",
   "metadata": {},
   "outputs": [],
   "source": [
    "def spec_rad(M):\n",
    "    \"\"\"\n",
    "    Compute the spectral radius of M.\n",
    "    \"\"\"\n",
    "    return np.max(np.abs(np.linalg.eigvals(M)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bdca499",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[1,2],[2,1]])\n",
    "spec_rad(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4fa90a5",
   "metadata": {},
   "source": [
    "This function is available in the `quantecon_book_networks` package, along with \n",
    "several other functions for used repeatedly in the text.  Source code for\n",
    "these functions can be seen [here](pkg_funcs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a1e4cba",
   "metadata": {},
   "outputs": [],
   "source": [
    "qbn_io.spec_rad(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ac85796",
   "metadata": {},
   "source": [
    "## Probability\n",
    "\n",
    "### The unit simplex in $\\mathbb{R}^3$\n",
    "\n",
    "Here we define a function for plotting the unit simplex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a661d596",
   "metadata": {},
   "outputs": [],
   "source": [
    "def unit_simplex(angle):\n",
    "    \n",
    "    fig = plt.figure(figsize=(10, 8))\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "    vtx = [[0, 0, 1],\n",
    "           [0, 1, 0], \n",
    "           [1, 0, 0]]\n",
    "    \n",
    "    tri = Poly3DCollection([vtx], color='darkblue', alpha=0.3)\n",
    "    tri.set_facecolor([0.5, 0.5, 1])\n",
    "    ax.add_collection3d(tri)\n",
    "\n",
    "    ax.set(xlim=(0, 1), ylim=(0, 1), zlim=(0, 1), \n",
    "           xticks=(1,), yticks=(1,), zticks=(1,))\n",
    "\n",
    "    ax.set_xticklabels(['$(1, 0, 0)$'], fontsize=16)\n",
    "    ax.set_yticklabels([f'$(0, 1, 0)$'], fontsize=16)\n",
    "    ax.set_zticklabels([f'$(0, 0, 1)$'], fontsize=16)\n",
    "\n",
    "    ax.xaxis.majorTicks[0].set_pad(15)\n",
    "    ax.yaxis.majorTicks[0].set_pad(15)\n",
    "    ax.zaxis.majorTicks[0].set_pad(35)\n",
    "\n",
    "    ax.view_init(30, angle)\n",
    "\n",
    "    # Move axis to origin\n",
    "    ax.xaxis._axinfo['juggled'] = (0, 0, 0)\n",
    "    ax.yaxis._axinfo['juggled'] = (1, 1, 1)\n",
    "    ax.zaxis._axinfo['juggled'] = (2, 2, 0)\n",
    "    \n",
    "    ax.grid(False) \n",
    "    \n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33c4089d",
   "metadata": {},
   "source": [
    "We can now produce the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "337d3792",
   "metadata": {},
   "outputs": [],
   "source": [
    "unit_simplex(50)\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/simplex_1.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8086d97f",
   "metadata": {},
   "source": [
    "### Independent draws from Student’s t and Normal distributions\n",
    "\n",
    "Here we illustrate the occurrence of \"extreme\" events in heavy tailed distributions. \n",
    "We start by generating 1,000 samples from a normal distribution and a Student's t distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7be0e011",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import t\n",
    "n = 1000\n",
    "np.random.seed(123)\n",
    "\n",
    "s = 2\n",
    "n_data = np.random.randn(n) * s\n",
    "\n",
    "t_dist = t(df=1.5)\n",
    "t_data = t_dist.rvs(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37b4b6cb",
   "metadata": {},
   "source": [
    "When we plot our samples, we see the Student's t distribution frequently\n",
    "generates samples many standard deviations from the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90cd4891",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(8, 3.4))\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_ylim((-50, 50))\n",
    "    ax.plot((0, n), (0, 0), 'k-', lw=0.3)\n",
    "\n",
    "ax = axes[0]\n",
    "ax.plot(list(range(n)), t_data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(list(range(n)), 0, t_data, 'k', lw=0.2)\n",
    "ax.get_xaxis().set_major_formatter(\n",
    "    ticker.FuncFormatter(lambda x, p: format(int(x), ',')))\n",
    "ax.set_title(f\"Student t draws\", fontsize=11)\n",
    "\n",
    "ax = axes[1]\n",
    "ax.plot(list(range(n)), n_data, linestyle='', marker='o', alpha=0.5, ms=4)\n",
    "ax.vlines(list(range(n)), 0, n_data, lw=0.2)\n",
    "ax.get_xaxis().set_major_formatter(\n",
    "    ticker.FuncFormatter(lambda x, p: format(int(x), ',')))\n",
    "ax.set_title(f\"$N(0, \\sigma^2)$ with $\\sigma = {s}$\", fontsize=11)\n",
    "\n",
    "plt.tight_layout()\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/heavy_tailed_draws.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca8d8028",
   "metadata": {},
   "source": [
    "### CCDF plots for the Pareto and Exponential distributions\n",
    "\n",
    "When the Pareto tail property holds, the CCDF is eventually log linear. Here\n",
    "we illustrates this using a Pareto distribution. For comparison, an exponential\n",
    "distribution is also shown. First we define our domain and the Pareto and\n",
    "Exponential distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a79084d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(1, 10, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89d0baa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "α = 1.5\n",
    "def Gp(x):\n",
    "    return x**(-α)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb412bf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "λ = 1.0\n",
    "def Ge(x):\n",
    "    return np.exp(-λ * x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90785fac",
   "metadata": {},
   "source": [
    "We then plot our distribution on a log-log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6b2f844",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=default_figsize)\n",
    "\n",
    "ax.plot(np.log(x), np.log(Gp(x)), label=\"Pareto\")\n",
    "ax.plot(np.log(x), np.log(Ge(x)), label=\"Exponential\")\n",
    "\n",
    "ax.legend(fontsize=12, frameon=False, loc=\"lower left\")\n",
    "ax.set_xlabel(\"$\\ln x$\", fontsize=12)\n",
    "ax.set_ylabel(\"$\\ln G(x)$\", fontsize=12)\n",
    "\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/ccdf_comparison_1.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63906e68",
   "metadata": {},
   "source": [
    "### Empirical CCDF plots for largest firms (Forbes)\n",
    "\n",
    "Here we show that the distribution of firm sizes has a Pareto tail. We start\n",
    "by loading the `forbes_global_2000` dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b9e33a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfff = ch1_data['forbes_global_2000']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9266e970",
   "metadata": {},
   "source": [
    "We calculate values of the empirical CCDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a342d5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.asarray(dfff['Market Value'])[0:500]\n",
    "y_vals = np.empty_like(data, dtype='float64')\n",
    "n = len(data)\n",
    "for i, d in enumerate(data):\n",
    "    # record fraction of sample above d\n",
    "    y_vals[i] = np.sum(data >= d) / n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "690fa754",
   "metadata": {},
   "source": [
    "Now we fit a linear trend line (on the log-log scale)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "956c019d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.log(data), np.log(y_vals)\n",
    "results = sm.OLS(y, sm.add_constant(x)).fit()\n",
    "b, a = results.params"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c90d836b",
   "metadata": {},
   "source": [
    "Finally we produce our plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a13d4b05",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(7.3, 4))\n",
    "\n",
    "ax.scatter(x, y, alpha=0.3, label=\"firm size (market value)\")\n",
    "ax.plot(x, x * a + b, 'k-', alpha=0.6, label=f\"slope = ${a: 1.2f}$\")\n",
    "\n",
    "ax.set_xlabel('log value', fontsize=12)\n",
    "ax.set_ylabel(\"log prob.\", fontsize=12)\n",
    "ax.legend(loc='lower left', fontsize=12)\n",
    "    \n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/empirical_powerlaw_plots_firms_forbes.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d80c5a0a",
   "metadata": {},
   "source": [
    "## Graph Theory\n",
    "\n",
    "### Zeta and Pareto distributions\n",
    "\n",
    "We begin by defining the Zeta and Pareto distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19032324",
   "metadata": {},
   "outputs": [],
   "source": [
    "γ = 2.0\n",
    "α = γ - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abccb3b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def z(k, c=2.0):\n",
    "    return c * k**(-γ)\n",
    "\n",
    "k_grid = np.arange(1, 10+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb68c8bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def p(x, c=2.0):\n",
    "    return c * x**(-γ)\n",
    "\n",
    "x_grid = np.linspace(1, 10, 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f605c4dc",
   "metadata": {},
   "source": [
    "Then we can produce our plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3661175",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=default_figsize)\n",
    "ax.plot(k_grid, z(k_grid), '-o', label='zeta distribution with $\\gamma=2$')\n",
    "ax.plot(x_grid, p(x_grid), label='density of Pareto with tail index $\\\\alpha$')\n",
    "ax.legend(fontsize=12)\n",
    "ax.set_yticks((0, 1, 2))\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/zeta_1.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35f01cce",
   "metadata": {},
   "source": [
    "### NetworkX digraph plot\n",
    "\n",
    "We start by creating a graph object and populating it with edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd8fae62",
   "metadata": {},
   "outputs": [],
   "source": [
    "G_p = nx.DiGraph()\n",
    "\n",
    "edge_list = [\n",
    "    ('p', 'p'),\n",
    "    ('m', 'p'), ('m', 'm'), ('m', 'r'),\n",
    "    ('r', 'p'), ('r', 'm'), ('r', 'r')\n",
    "]\n",
    "\n",
    "for e in edge_list:\n",
    "    u, v = e\n",
    "    G_p.add_edge(u, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77e32349",
   "metadata": {},
   "source": [
    "Now we can plot our graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "892bc151",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "nx.spring_layout(G_p, seed=4)\n",
    "nx.draw_spring(G_p, ax=ax, node_size=500, with_labels=True, \n",
    "                 font_weight='bold', arrows=True, alpha=0.8,\n",
    "                 connectionstyle='arc3,rad=0.25', arrowsize=20)\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/networkx_basics_1.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f089c248",
   "metadata": {},
   "source": [
    "The `DiGraph` object has methods that calculate in-degree and out-degree of vertices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3e1adc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "G_p.in_degree('p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a502c57a",
   "metadata": {},
   "outputs": [],
   "source": [
    "G_p.out_degree('p')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84636840",
   "metadata": {},
   "source": [
    "Additionally, the `NetworkX` package supplies functions for testing\n",
    "communication and strong connectedness, as well as to compute strongly\n",
    "connected components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b41f9be",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.DiGraph()\n",
    "G.add_edge(1, 1)\n",
    "G.add_edge(2, 1)\n",
    "G.add_edge(2, 3)\n",
    "G.add_edge(3, 2)\n",
    "list(nx.strongly_connected_components(G))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79d4b29c",
   "metadata": {},
   "source": [
    "Like `NetworkX`, the Python library `quantecon` \n",
    "provides access to some graph-theoretic algorithms. \n",
    "\n",
    "In the case of QuantEcon's `DiGraph` object, an instance is created via the adjacency matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbd43d1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = ((1, 0, 0),\n",
    "     (1, 1, 1),\n",
    "     (1, 1, 1))\n",
    "A = np.array(A) # Convert to NumPy array\n",
    "G = qe.DiGraph(A)\n",
    "\n",
    "G.strongly_connected_components"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c670e09",
   "metadata": {},
   "source": [
    "### International private credit flows by country\n",
    "\n",
    "We begin by loading an adjacency matrix of international private credit flows\n",
    "(in the form of a NumPy array and a list of country labels)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd6f046d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Z = ch1_data[\"adjacency_matrix\"][\"Z\"]\n",
    "Z_visual= ch1_data[\"adjacency_matrix\"][\"Z_visual\"]\n",
    "countries = ch1_data[\"adjacency_matrix\"][\"countries\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ade862b8",
   "metadata": {},
   "source": [
    "To calculate our graph's properties, we use hub-based eigenvector\n",
    "centrality as our centrality measure for this plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a89769e",
   "metadata": {},
   "outputs": [],
   "source": [
    "centrality = qbn_io.eigenvector_centrality(Z_visual, authority=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7be043ae",
   "metadata": {},
   "source": [
    "Now we convert our graph features to plot features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce2696fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "node_colors = cm.plasma(qbn_io.to_zero_one_beta(centrality))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ba53d05",
   "metadata": {},
   "source": [
    "Finally we produce the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f802d18",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = qbn_io.to_zero_one_beta(Z.sum(axis=1))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 10))\n",
    "plt.axis(\"off\")\n",
    "\n",
    "qbn_plot.plot_graph(Z_visual, X, ax, countries,\n",
    "           layout_type='spring',\n",
    "           layout_seed=1234,\n",
    "           node_size_multiple=3000,\n",
    "           edge_size_multiple=0.000006,\n",
    "           tol=0.0,\n",
    "           node_color_list=node_colors) \n",
    "\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/financial_network_analysis_visualization.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88e9c96e",
   "metadata": {},
   "source": [
    "### Centrality measures for the credit network\n",
    "\n",
    "This figure looks at six different centrality measures.\n",
    "\n",
    "We begin by defining a function for calculating eigenvector centrality.\n",
    "\n",
    "Hub-based centrality is calculated by default, although authority-based centrality\n",
    "can be calculated by setting `authority=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "577372e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def eigenvector_centrality(A, k=40, authority=False):\n",
    "    \"\"\"\n",
    "    Computes the dominant eigenvector of A. Assumes A is \n",
    "    primitive and uses the power method.  \n",
    "    \n",
    "    \"\"\"\n",
    "    A_temp = A.T if authority else A\n",
    "    n = len(A_temp)\n",
    "    r = spec_rad(A_temp)\n",
    "    e = r**(-k) * (np.linalg.matrix_power(A_temp, k) @ np.ones(n))\n",
    "    return e / np.sum(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16682dc",
   "metadata": {},
   "source": [
    "Here a similar function is defined for calculating Katz centrality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d7fbee4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def katz_centrality(A, b=1, authority=False):\n",
    "    \"\"\"\n",
    "    Computes the Katz centrality of A, defined as the x solving\n",
    "\n",
    "    x = 1 + b A x    (1 = vector of ones)\n",
    "\n",
    "    Assumes that A is square.\n",
    "\n",
    "    If authority=True, then A is replaced by its transpose.\n",
    "    \"\"\"\n",
    "    n = len(A)\n",
    "    I = np.identity(n)\n",
    "    C = I - b * A.T if authority else I - b * A\n",
    "    return np.linalg.solve(C, np.ones(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0812961",
   "metadata": {},
   "source": [
    "Now we generate an unweighted version of our matrix to help calculate in-degree and out-degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d665a8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "D = qbn_io.build_unweighted_matrix(Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "947d8b01",
   "metadata": {},
   "source": [
    "We now use the above to calculate the six centrality measures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "413f7029",
   "metadata": {},
   "outputs": [],
   "source": [
    "outdegree = D.sum(axis=1)\n",
    "ecentral_hub = eigenvector_centrality(Z, authority=False)\n",
    "kcentral_hub = katz_centrality(Z, b=1/1_700_000)\n",
    "\n",
    "indegree = D.sum(axis=0)\n",
    "ecentral_authority = eigenvector_centrality(Z, authority=True)\n",
    "kcentral_authority = katz_centrality(Z, b=1/1_700_000, authority=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a14bf02e",
   "metadata": {},
   "source": [
    "Here we provide a helper function that returns a DataFrame for each measure.\n",
    "The DataFrame is ordered by that measure and contains color information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bcb86db",
   "metadata": {},
   "outputs": [],
   "source": [
    "def centrality_plot_data(countries, centrality_measures):\n",
    "    df = pd.DataFrame({'code': countries,\n",
    "                       'centrality':centrality_measures, \n",
    "                       'color': qbn_io.colorise_weights(centrality_measures).tolist()\n",
    "                       })\n",
    "    return df.sort_values('centrality')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5777a24",
   "metadata": {},
   "source": [
    "Finally, we plot the various centrality measures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47130789",
   "metadata": {},
   "outputs": [],
   "source": [
    "centrality_measures = [outdegree, indegree, \n",
    "                       ecentral_hub, ecentral_authority, \n",
    "                       kcentral_hub, kcentral_authority]\n",
    "\n",
    "ylabels = ['out degree', 'in degree',\n",
    "           'eigenvector hub','eigenvector authority', \n",
    "           'Katz hub', 'Katz authority']\n",
    "\n",
    "ylims = [(0, 20), (0, 20), \n",
    "         None, None,   \n",
    "         None, None]\n",
    "\n",
    "\n",
    "fig, axes = plt.subplots(3, 2, figsize=(10, 12))\n",
    "\n",
    "axes = axes.flatten()\n",
    "\n",
    "for i, ax in enumerate(axes):\n",
    "    df = centrality_plot_data(countries, centrality_measures[i])\n",
    "      \n",
    "    ax.bar('code', 'centrality', data=df, color=df[\"color\"], alpha=0.6)\n",
    "    \n",
    "    patch = mpatches.Patch(color=None, label=ylabels[i], visible=False)\n",
    "    ax.legend(handles=[patch], fontsize=12, loc=\"upper left\", handlelength=0, frameon=False)\n",
    "    \n",
    "    if ylims[i] is not None:\n",
    "        ax.set_ylim(ylims[i])\n",
    "\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/financial_network_analysis_centrality.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11204e33",
   "metadata": {},
   "source": [
    "### Computing in and out degree distributions\n",
    "\n",
    "The in-degree distribution evaluated at $k$ is the fraction of nodes in a\n",
    "network that have in-degree $k$. The in-degree distribution of a `NetworkX`\n",
    "DiGraph can be calculated using the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61512e1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def in_degree_dist(G):\n",
    "    n = G.number_of_nodes()\n",
    "    iG = np.array([G.in_degree(v) for v in G.nodes()])\n",
    "    d = [np.mean(iG == k) for k in range(n+1)]\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8121ffff",
   "metadata": {},
   "source": [
    "The out-degree distribution is defined analogously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc3e776b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def out_degree_dist(G):\n",
    "    n = G.number_of_nodes()\n",
    "    oG = np.array([G.out_degree(v) for v in G.nodes()])\n",
    "    d = [np.mean(oG == k) for k in range(n+1)]\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "840131af",
   "metadata": {},
   "source": [
    "### Degree distribution for international aircraft trade\n",
    "\n",
    "Here we illustrate that the commercial aircraft international trade network is\n",
    "approximately scale-free by plotting the degree distribution alongside\n",
    "$f(x)=cx-\\gamma$ with $c=0.2$ and $\\gamma=1.1$. \n",
    "\n",
    "In this calculation of the degree distribution, performed by the NetworkX\n",
    "function `degree_histogram`, directions are ignored and the network is treated\n",
    "as an undirected graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecc6928e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_degree_dist(G, ax, loglog=True, label=None):\n",
    "    \"Plot the degree distribution of a graph G on axis ax.\"\n",
    "    dd = [x for x in nx.degree_histogram(G) if x > 0]\n",
    "    dd = np.array(dd) / np.sum(dd)  # normalize\n",
    "    if loglog:\n",
    "        ax.loglog(dd, '-o', lw=0.5, label=label)\n",
    "    else:\n",
    "        ax.plot(dd, '-o', lw=0.5, label=label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5ec0060",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=default_figsize)\n",
    "\n",
    "plot_degree_dist(DG, ax, loglog=False, label='degree distribution')\n",
    "\n",
    "xg = np.linspace(0.5, 25, 250)\n",
    "ax.plot(xg, 0.2 * xg**(-1.1), label='power law')\n",
    "ax.set_xlim(0.9, 22)\n",
    "ax.set_ylim(0, 0.25)\n",
    "ax.legend()\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/commercial_aircraft_2019_2.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c8cb9ba",
   "metadata": {},
   "source": [
    "### Random graphs\n",
    "\n",
    "The code to produce the Erdos-Renyi random graph, used below, applies the\n",
    "combinations function from the `itertools` library. The function\n",
    "`combinations(A, k)` returns a list of all subsets of $A$ of size $k$. For\n",
    "example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "875a6ae1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "letters = 'a', 'b', 'c'\n",
    "list(itertools.combinations(letters, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe83b981",
   "metadata": {},
   "source": [
    "Below we generate random graphs using the Erdos-Renyi and Barabasi-Albert\n",
    "algorithms. Here, for convenience, we will define a function to plot these\n",
    "graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc9858c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_random_graph(RG,ax):\n",
    "    node_pos_dict = nx.spring_layout(RG, k=1.1)\n",
    "\n",
    "    centrality = nx.degree_centrality(RG)\n",
    "    node_color_list = qbn_io.colorise_weights(list(centrality.values()))\n",
    "\n",
    "    edge_color_list = []\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            edge_color_list.append(node_color_list[i])\n",
    "\n",
    "    nx.draw_networkx_nodes(RG, \n",
    "                           node_pos_dict, \n",
    "                           node_color=node_color_list, \n",
    "                           edgecolors='grey', \n",
    "                           node_size=100,\n",
    "                           linewidths=2, \n",
    "                           alpha=0.8, \n",
    "                           ax=ax)\n",
    "\n",
    "    nx.draw_networkx_edges(RG, \n",
    "                           node_pos_dict, \n",
    "                           edge_color=edge_colors, \n",
    "                           alpha=0.4,  \n",
    "                           ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "968c2f98",
   "metadata": {},
   "source": [
    "### An instance of an Erdos–Renyi random graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "027cd2dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "p = 0.05\n",
    "G_er = qbn_io.erdos_renyi_graph(n, p, seed=1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86c89612",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.2))\n",
    "\n",
    "axes[0].set_title(\"Graph visualization\")\n",
    "plot_random_graph(G_er,axes[0])\n",
    "\n",
    "axes[1].set_title(\"Degree distribution\")\n",
    "plot_degree_dist(G_er, axes[1], loglog=False)\n",
    "\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/rand_graph_experiments_1.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c2f0411",
   "metadata": {},
   "source": [
    "### An instance of a preferential attachment random graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fcae0d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "m = 5\n",
    "G_ba = nx.generators.random_graphs.barabasi_albert_graph(n, m, seed=123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b513178",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.2))\n",
    "\n",
    "axes[0].set_title(\"Graph visualization\")\n",
    "plot_random_graph(G_ba, axes[0])\n",
    "\n",
    "axes[1].set_title(\"Degree distribution\")\n",
    "plot_degree_dist(G_ba, axes[1], loglog=False)\n",
    "\n",
    "if export_figures:\n",
    "    plt.savefig(\"figures/rand_graph_experiments_2.pdf\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
